% X \in S^{n(p_max+1)}
% C \in S^{n(p_max+1)}
% TODO decompose into functions? not sure how to do that without losing cvx scope
% TODO vectorize; this works but it slow
function B = regular_mlm(C, n, p_max, gamma, lambda_0)
	cvx_begin
		variable X(n*(p_max+1), n*(p_max+1)) symmetric

% construct D
		expression D(n, n*(p_max+1));
		for i=0:p_max
			D(1:n,1:n) = D(1:n,1:n) + X((1:n)+n*i,(1:n)+n*i);
        end		
        for k=1:p_max
			for i=0:(p_max-k)
				D(1:n,(1:n)+n*k) = D(1:n,(1:n)+n*k) + 2*X((1:n)+n*i,(1:n)+n*(i+k));
			end
        end

% evaluate penalty h
		expression h(1);
		for i=1:n
			for j=(i+1):n
				h = h + norm([D(i,j+n*(0:p_max)) D(j,i+n*(0:p_max))], Inf);
			end
		end

% evaluate block penalty g
% TODO try different penalties g
		expression g(1);
		for i=1:p_max
% penalty parameter is linear in index of block; TODO try logarithmic?
%			lambda_i = lambda_0 * i;
%			lambda_i = lambda_0 * log(i);
			lambda_i = lambda_0;
			for j=(i+1):p_max
				g = g + lambda_i * max(max(abs(X((1:n)+n*i,(1:n)+n*j))));
			end
		end

		minimize(-log_det(X(1:n,1:n)) + trace(C*X) + gamma*h + g);
		X == semidefinite(n*(p_max+1));
	cvx_end

	B0 = sqrtm(X(1:n,1:n));
	B1p = B0 \ X(1:n,(1+n):(n*(p_max+1)));
	B = [B0 B1p];
end
